addpath('/Users/nicolegorton/Dropbox/World Bank/OTN_LAC/proj_code/Toolbox/table2latex/');

folders;

% This file loads in baseline grids and OSM grids, compares welfare gains
% and location of investments

country_list_short;
CFAC_TYPE = {'exp_eng', 'mis_eng'}; %{'exp_eng','exp_til','mis_eng'};

% parameters for counterfactual grid
adj_width = 6;
max_bright = 1.5;
gamma = 0.10;
mobil = 'off'; % 'on',
cong = 'on';
alpha = 0.4;
beta = 0.13;
sigma = 5;
rho = 0;
a = 1;
Ngoods = 10;
nu = 1;

ROBUST = {'BASE', 'WP', 'OSM', 'WP-OSM'};

param = init_parameters( 'a',a,'rho',rho,'alpha',alpha,'sigma',sigma,...       % preferences and technology
    'beta',beta,'gamma',gamma,'nu',nu,'m',ones(Ngoods+1,1),...   % transpot costs
    'K',1,...
    'LaborMobility',char(mobil),...
    'N',Ngoods+1,...
    'CrossGoodCongestion',char(cong),...
    'TolKappa',1e-4 );

% Initiate a blank table for results
%{ 
welfare_exp_pre = zeros(Ncountries, 1);
welfare_exp_post = zeros(Ncountries, 1);
welfare_exp_WP = zeros(Ncountries, 1);
welfare_mis_pre = zeros(Ncountries, 1);
welfare_mis_post = zeros(Ncountries, 1);
welfare_mis_WP = zeros(Ncountries, 1);
%}

welfare_exp = zeros(Ncountries, 4);
welfare_mis = zeros(Ncountries, 4);

corr_exp = zeros(Ncountries, 3);
corr_mis = zeros(Ncountries, 3);

results = table(country_names(:, 3), welfare_exp, welfare_mis, corr_exp, corr_mis);

for country_n = 1:Ncountries
    
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );
    
    % Loop through counterfactuals: misallocation and expansion
    for cfac_type = CFAC_TYPE        

        default_cellsize = 0.5;                    

        x_ver_hor = get_threshold( country ) ;   
        x_diag = x_ver_hor;                        
        
        % load in the country graph (shouldn't matter which one really)
        graph_filename = [ path_load_grids,country,'_grid_',...
        num2str( x_diag ),'_',...
        num2str( x_ver_hor ), '_', num2str(default_cellsize), ...
        '_connected1_LAC_BASE.mat' ];

        load( graph_filename )
        unique_edges=country_graph.unique_edges;
            
        rcount = 1;
        for robust = ROBUST 
             
            [ robust ]
            if strcmp(robust, 'BASE')
                pop_source = 0;
                osm = 0;      
            elseif strcmp(robust, 'WP')
                pop_source = 1;
                osm = 0;
            elseif strcmp(robust, 'OSM')
                pop_source = 0;
                osm = 1;          
            elseif strcmp(robust, 'WP-OSM')
                pop_source = 1;
                osm = 1;             
            end   

            % define file to load calibration and set diary
            filename = [ country,...
                        '_diag',num2str( x_diag ),...
                        '_hor',num2str( x_ver_hor ),...
                        '_a',num2str( param.a ),...
                        '_rho',num2str( param.rho ),...
                        '_alpha',num2str( param.alpha ),...
                        '_sigma',num2str( param.sigma ),...
                        '_beta',num2str( param.beta ),...
                        '_gamma',num2str( param.gamma ),...
                        '_nu',num2str( param.nu ),...
                        '_mobil',num2str( param.mobility ),...
                        '_cong',num2str( param.cong ),...
                        '_ngoods',num2str( Ngoods ),...
                        '_cellsize',num2str(default_cellsize) ,...
                        '_osm', num2str(osm), '_WP', num2str(pop_source)];

        
         if  ~strcmp(country, 'Brazil') || (strcmp(country, 'Brazil') && ismember(robust, {'BASE', 'OSM'})) 
             
            % load the counterfactual grid 
            load( [ path_save_cfactuals,filename,'_cfactual_', char(cfac_type), '_' , char(robust) , '.mat' ],'cfactual' );                       
            
            [cfactual.welfare_gain]
            welfare = (cfactual.welfare_gain-1)*100;
            if strcmp(string(cfac_type), 'exp_eng')
                if  welfare<1 
                    results.welfare_exp(country_n ,rcount) = round(welfare, 1);
                else
                    results.welfare_exp(country_n ,rcount) = round(welfare, 1);                                       
                end
            elseif strcmp(string(cfac_type), 'mis_eng')
                if  welfare<1 
                    results.welfare_mis(country_n ,rcount) = round(welfare, 1);
                else
                    results.welfare_mis(country_n ,rcount) = round(welfare, 1);                                       
               end
            end
            
            if strcmp(robust, 'BASE')
                base_cfactual=cfactual;
                I_obs = zeros( length(unique_edges),1 );
                I_exp = zeros( length(unique_edges),1 );

                for j=1:length(unique_edges)

                    I_obs(j) = base_cfactual.results_actual.Ijk( unique_edges(j,1),unique_edges(j,2) );
                    I_exp(j) = base_cfactual.results_cfactual.Ijk( unique_edges(j,1),unique_edges(j,2) );

                end    
                dI = I_exp-I_obs;                                            
                
            else
                
                % lanes in calibrated and counterfactual model
                I_obs_new = zeros( length(unique_edges),1 );
                I_exp_new = zeros( length(unique_edges),1 );
                for j=1:length(unique_edges)

                    I_obs_new(j) = cfactual.results_actual.Ijk( unique_edges(j,1),unique_edges(j,2) );
                    I_exp_new(j) = cfactual.results_cfactual.Ijk( unique_edges(j,1),unique_edges(j,2) );

                end

                % compute changes
                dI_new = I_exp_new-I_obs_new;                           
            
                if strcmp(string(cfac_type), 'exp_eng')
                    results.corr_exp(country_n, rcount-1) = round(corr(dI, dI_new)*100, 1);
                else
                    results.corr_mis(country_n, rcount-1) = round(corr(dI, dI_new)*100, 1);
                end
            end
             
         end
            rcount = rcount+1;            
        end
    end
end

results = sortrows(results,2, 'descend')
writetable(results,  [path_notes_images, 'robustness.xlsx'])

